1789F - Serval and Brain Power - CodeForces Solution

#include <bits/stdc++.h>

#pragma region
namespace zayin {

namespace io {
struct convertor {
  template <typename, typename = void>
  struct is_iterable {
    static constexpr bool value = false;
  };
  template <typename T>
  struct is_iterable<
      T,
      std::void_t<decltype(std::declval<T>().begin()), decltype(std::declval<T>().end())>> {
    static constexpr bool value = true;
  };

  template <typename, typename = void>
  struct is_to_string_implemented {
    static constexpr bool value = false;
  };
  template <typename T>
  struct is_to_string_implemented<T, std::void_t<decltype(std::declval<T>().to_string())>> {
    static constexpr bool value = true;
  };

  template <typename F, typename S>
  static std::string to_string(const std::pair<F, S>& p) {
    return "(" + to_string(p.first) + ", " + to_string(p.second) + ")";
  }
  static std::string to_string(const std::string& s) { return '"' + s + '"'; }
  static std::string to_string(const char* s) { return '"' + std::string(s) + '"'; }
  static std::string to_string(const char c) { return "'" + std::string(1, c) + "'"; }
  static std::string to_string(const bool b) { return std::string(1, "FT"[b]); }
  template <typename AT, typename std::enable_if_t<std::is_arithmetic<AT>::value, int> = 0>
  static std::string to_string(const AT a) {
    return std::to_string(a);
  }

  template <typename VT, typename std::enable_if_t<is_iterable<VT>::value, int> = 0>
  static std::string to_string(const VT& v, const int indent_width = 0) {
    if (!v.size()) return "{}";
    std::string res;
    if constexpr (is_iterable<decltype(*v.begin())>::value) {
      res = to_string(*v.begin(), indent_width + 1);
      for (auto it = std::next(v.begin()); it != v.end(); ++it)
        res += ",\n" + to_string(*it, indent_width + 1);
      res[indent_width] = '{';
    } else {
      res = std::string(indent_width, ' ') + "{" + to_string(*v.begin());
      for (auto it = std::next(v.begin()); it != v.end(); ++it) res += ", " + to_string(*it);
    }
    res += "}";
    return res;
  }

  template <typename T, typename std::enable_if_t<is_to_string_implemented<T>::value, int> = 0>
  static std::string to_string(T t) {
    return t.to_string();
  }
};

class ostream {
 public:
  ostream(FILE* file) : file(file) {}
  ~ostream() { reflesh(); }

  template <typename T>
  ostream& operator<<(T t) {
    write_single(t);
    return *this;
  }

  template <char... escapes>
  void log() {
    write_single('\n');
  }
  template <char... escapes, typename T, typename... Args>
  void log(T t, Args... args) {
    write_single(t);
    if (sizeof...(args)) (write_single(escapes), ...);
    log<escapes...>(args...);
  }
  template <typename... Args>
  void operator()(Args... args) {
    log<' '>(args...);
  }
  void flush() {
    reflesh();
    fflush(file);
  }
  void reflesh() {
    fwrite(buffer, sizeof(char), pointer - buffer, file);
    pointer = buffer;
  }

 private:
  void write_single(const char& c) {
    if (pointer == buffer + buffer_size) reflesh();
    *(pointer++) = c;
  }
  void write_single(const char* s) { write_single(std::string(s)); }
  void write_single(const std::string& s) {
    for (char c : s) write_single(c);
  }
  template <typename T>
  void write_single(T t) {
    for (char c : convertor::to_string(t)) write_single(c);
  }
  FILE* file;
  static const int buffer_size = 1 << 23;
  char buffer[buffer_size], *pointer = buffer;
};
ostream cout(stdout);
ostream cerr(stderr);
const char endl = '\n';

class istream : public std::istream {
 public:
  istream(std::istream& is) : is(is) {}
  template <typename T>
  istream& operator>>(T& t) {
    is >> t;
    return *this;
  }
  void operator()() {}
  template <typename T, typename... Args>
  void operator()(T& t, Args&... args) {
    is >> t;
    operator()(args...);
  }

  operator bool() { return (bool)is; }

 private:
  std::istream& is;
};
istream cin(std::cin);

class debug_helper {
 public:
  static bool isLeftBracket(char c) { return c == '[' || c == '(' || c == '{'; }
  static bool isRightBracket(char c) { return c == ']' || c == ')' || c == '}'; }
  template <typename ostreamType, typename... Args>
  static void dump(ostreamType& os, std::string names, Args... args) {
    std::vector<std::string> val;
    (val.push_back(convertor::to_string(args)), ...);

    std::vector<std::string> name{""};
    int count = 0;
    for (char c : names) {
      if (isspace(c)) continue;
      if (c == ',' && !count) {
        name.push_back("");
        continue;
      }
      if (isLeftBracket(c)) ++count;
      if (isRightBracket(c)) --count;
      name.back() += c;
    }
    assert(name.size() == val.size());
    std::unordered_map<std::string, std::string> mappings;
    for (int i = 0; i < name.size(); ++i) mappings[name[i]] = val[i];
    auto get_val = [&](std::string name) -> std::optional<std::string> {
      if (mappings.count(name)) return mappings[name];
      for (auto pr : mappings) {
        for (int i = 0; i + pr.first.size() <= name.size(); ++i) {
          if (name.substr(i, pr.first.size()) != pr.first) continue;
          if (i && std::isalnum(name[i - 1])) continue;
          if (i + 1 < name.size() && std::isalnum(name[i + 1])) continue;
          name.replace(i, pr.first.size(), pr.second);
        }
      }
      return get_expression_val(name);
    };
    std::function<std::string(std::string)> simplify = [&](std::string name) -> std::string {
      assert(mappings.count(name));
      std::stack<std::string> stk;
      stk.push("");
      bool should_add_mapping = true;
      for (char c : name) {
        if (isLeftBracket(c)) {
          stk.top() += c, stk.push("");
        } else if (isRightBracket(c) || c == ',') {
          std::optional<std::string> val = get_val(stk.top());
          if (!val.has_value()) {
            should_add_mapping = false;
            constexpr int kDisplayLength = 5;
            val = stk.top().size() > kDisplayLength ? "*" : stk.top();
          }
          stk.pop();
          stk.top() += val.value() + c;
        } else {
          stk.top() += c;
        }
      }
      assert(stk.size() == 1);
      if (should_add_mapping) mappings[stk.top()] = mappings[name];
      return stk.top();
    };
    for (int i = 0; i < name.size(); ++i) {
      if (i) os << ", ";
      os << simplify(name[i]) << " = " << mappings.at(name[i]);
    }
    os << '\n';
    os.reflesh();
  }

 private:
  static std::optional<std::string> get_expression_val(std::string expression) {
    if (expression.empty()) return "";
    if (std::any_of(expression.begin(), expression.end(), [&](char c) {
          if ('0' <= c && c <= '9') return false;
          if (c == '+' || c == '-' || c == '*' || c == '(' || c == ')') return false;
          return true;
        }))
      return std::nullopt;
    std::stack<long long> num_stk;
    std::stack<long long> op_stk;
    auto get_precedence = [&](char c) {
      if (c == '#') return std::numeric_limits<int>::min();
      if (c == '(') return -2;
      if (c == ')') return -1;
      if (c == '+' || c == '-') return 0;
      if (c == '*') return 1;
    };
    auto get_op_result = [&](long long a, char op, long long b) {
      if (op == '+') return a + b;
      if (op == '-') return a - b;
      if (op == '*') return a * b;
    };
    expression += '#';
    for (int i = 0; i < expression.size(); ++i) {
      if (op_stk.size() == num_stk.size()) {
        long long sign = 1;
        if (!std::isalnum(expression[i])) {
          if (expression[i] == '+')
            sign = 1;
          else if (expression[i] == '-')
            sign = -1;
          else
            assert(0);
          ++i;
        }
        assert(std::isalnum(expression[i]));
        num_stk.push(0);
        for (; i < expression.size() && std::isalnum(expression[i]); ++i)
          num_stk.top() = num_stk.top() * 10 + expression[i] - '0';
        num_stk.top() = sign * num_stk.top();
        --i;
      } else if (expression[i] == '(') {
        op_stk.push('(');
      } else {
        while (op_stk.size() && get_precedence(op_stk.top()) >= get_precedence(expression[i])) {
          long long b = num_stk.top();
          num_stk.pop();
          long long a = num_stk.top();
          num_stk.pop();
          long long res = get_op_result(a, op_stk.top(), b);
          op_stk.pop();
          num_stk.push(res);
        }
        if (expression[i] == '+' || expression[i] == '-' || expression[i] == '*' ||
            expression[i] == '#') {
          op_stk.push(expression[i]);
        } else if (expression[i] == ')') {
          assert(op_stk.top() == '(');
          op_stk.pop();
        }
      }
    }
    assert(num_stk.size() == 1);
    assert(op_stk.size() == 1 && op_stk.top() == '#');
    return std::to_string(num_stk.top());
  };
};

#define dbg(args...) zayin::io::debug_helper::dump(zayin::io::cout, #args, ##args)

#define single_dump(x) #x << " = " << x
#define concat_dump(left, right) left << ", " << right
#define dump_to_stream(args...) dump_to_stream_id(args, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1)(args)
#define dump_to_stream_1(x1) single_dump(x1)
#define dump_to_stream_2(x1, x2) concat_dump(single_dump(x1), dump_to_stream_1(x2))
#define dump_to_stream_3(x1, x2, x3) concat_dump(single_dump(x1), dump_to_stream_2(x2, x3))
#define dump_to_stream_4(x1, x2, x3, x4) concat_dump(single_dump(x1), dump_to_stream_3(x2, x3, x4))
#define dump_to_stream_5(x1, x2, x3, x4, x5) \
  concat_dump(single_dump(x1), dump_to_stream_4(x2, x3, x4, x5))
#define dump_to_stream_6(x1, x2, x3, x4, x5, x6) \
  concat_dump(single_dump(x1), dump_to_stream_5(x2, x3, x4, x5, x6))
#define dump_to_stream_7(x1, x2, x3, x4, x5, x6, x7) \
  concat_dump(single_dump(x1), dump_to_stream_6(x2, x3, x4, x5, x6, x7))
#define dump_to_stream_8(x1, x2, x3, x4, x5, x6, x7, x8) \
  concat_dump(single_dump(x1), dump_to_stream_7(x2, x3, x4, x5, x6, x7, x8))
#define dump_to_stream_9(x1, x2, x3, x4, x5, x6, x7, x8, x9) \
  concat_dump(single_dump(x1), dump_to_stream_8(x2, x3, x4, x5, x6, x7, x8, x9))
#define dump_to_stream_10(x1, x2, x3, x4, x5, x6, x7, x8, x9, x10) \
  concat_dump(single_dump(x1), dump_to_stream_9(x2, x3, x4, x5, x6, x7, x8, x9, x10))
#define dump_to_stream_id(x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, args...) \
  dump_to_stream_##x11

}  // namespace io

namespace random {
static std::mt19937_64 engine(std::chrono::steady_clock::now().time_since_epoch().count());

template <typename T, typename std::enable_if_t<std::is_integral<T>::value, int> = 0>
T rand(T l, T r) {
  std::uniform_int_distribution<T> rng(l, r);
  return rng(engine);
}

}  // namespace random

namespace numerical {

long long slow_mul(long long a, long long b, long long modulus) {
  long long k = (long double)a * b / modulus;
  long long res = a * b - k * modulus;
  while (res < 0) res += modulus;
  while (res >= modulus) res -= modulus;
  return res;
}

namespace prime {

bool is_prime(long long x) {
  assert(x > 0);
  if (x == 1) return false;
  if (~x & 1) return x == 2;
  auto pow = [&](long long a, long long k) {
    long long res = 1;
    for (; k; k >>= 1) {
      if (k & 1) res = slow_mul(res, a, x);
      a = slow_mul(a, a, x);
    }
    return res;
  };
  int ctz = __builtin_ctzll(x - 1);
  auto miller_rabin = [&](long long a) {
    long long t = pow(a, (x - 1) >> ctz);
    for (int i = 0; i < ctz; ++i) {
      long long nt = slow_mul(t, t, x);
      if (nt == 1 && t != 1 && t != x - 1) return false;
      t = nt;
    }
    return t == 1;
  };
  static const std::vector<int> test_a = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
  for (int a : test_a) {
    if (x == a) return true;
    if (!miller_rabin(a)) return false;
  }
  return true;
}

// suppose x=prod pi^ei, return {pi,ei}
template <typename T, typename std::enable_if_t<std::is_integral<T>::value, int> = 0>
std::unordered_map<T, int> get_factorization(T x) {
  std::unordered_map<T, int> res;
  std::function<T(T)> pollard_rho = [&](T x) {
    while (1) {
      T c = random::rand<T>(0, x - 1);
      T x0 = random::rand<T>(1, x - 1), x1 = x0;
      T s = 1;
      for (int i = 1;; ++i) {
        x0 = (slow_mul(x0, x0, x) + c) % x;
        if (x0 == x1) break;
        s = slow_mul(s, std::abs(x0 - x1), x);
        T d = std::gcd(s, x);
        if (1 < d && d < x) return d;
        if (i > 1 && __builtin_popcount(i) == 1) x1 = x0;
      }
    }
  };
  std::function<void(T, int)> solve = [&](T x, int e) {
    if (x == 1) return;
    if (is_prime(x)) {
      res[x] += e;
      return;
    }
    T d = pollard_rho(x);
    int k = 0;
    for (; x % d == 0; x /= d) ++k;
    solve(x, e);
    solve(d, e * k);
  };
  solve(x, 1);
  return res;
}

template <typename T, typename std::enable_if_t<std::is_integral<T>::value, int> = 0>
std::vector<T> get_all_divisors(std::unordered_map<T, int> factorization) {
  std::vector<std::pair<T, int>> fv;
  for (auto [p, e] : factorization) fv.emplace_back(p, e);
  std::vector<T> res;
  std::function<void(int, T)> dfs = [&](int i, T d) {
    if (i == fv.size()) {
      res.push_back(d);
    } else {
      auto [p, max_e] = fv[i];
      for (int e = 0; e <= max_e; ++e) {
        dfs(i + 1, d);
        d *= p;
      }
    }
  };
  dfs(0, 1);
  std::sort(res.begin(), res.end());
  return res;
}

template <typename T, typename std::enable_if_t<std::is_integral<T>::value, int> = 0>
std::vector<T> get_all_divisors(T x) {
  return get_all_divisors(get_factorization(x));
}

}  // namespace prime

namespace modular {

struct modulus_type {
  static constexpr int root = 3;
  static constexpr int value = 998244353;
};

template <typename T>
class zint {
 public:
  using type = typename std::decay<decltype(T::value)>::type;

  constexpr static type modulus() { return T::value; }
  constexpr static zint root() { return zint(T::root); }

  constexpr zint() : value(0) {}
  template <typename U, typename std::enable_if_t<std::is_integral<U>::value, int> = 0>
  zint(const U& x) {
    value = normalize(x);
  }

  const type& operator()() const { return value; }
  template <typename U, typename std::enable_if_t<std::is_integral<U>::value, int> = 0>
  explicit operator U() const {
    return (U)value;
  }

  zint operator-() const { return zint(-value); }

  zint& operator+=(const zint& z) {
    value += z.value;
    if (value >= modulus()) value -= modulus();
    return *this;
  }
  zint& operator-=(const zint& z) {
    value -= z.value;
    if (value < 0) value += modulus();
    return *this;
  }
  zint& operator*=(const zint& z) {
    value = normalize((int64_t)value * (int64_t)z.value);
    return *this;
  }
  zint& operator/=(const zint& z) { return *this *= z.inv(); }
  zint pow(long long k) const {
    zint res = 1;
    for (zint z = *this; k; k >>= 1) {
      if (k & 1) res *= z;
      z *= z;
    }
    return res;
  }
  zint inv() const {
    assert(value != 0);
    type a = value, m = modulus();
    type u = 0, v = 1;
    while (a) {
      type t = m / a;
      m -= t * a;
      std::swap(a, m);
      u -= t * v;
      std::swap(u, v);
    }
    assert(m == 1);
    return zint(u);
  }
  zint operator+(const zint& z) const { return zint(*this) += z; }
  zint operator-(const zint& z) const { return zint(*this) -= z; }
  zint operator*(const zint& z) const { return zint(*this) *= z; }
  zint operator/(const zint& z) const { return zint(*this) /= z; }

  // TODO: support composite number
  bool is_sqrtable() const {
    if (!value) return true;
    return pow((modulus() - 1) / 2) == 1;
  }
  zint sqrt() const {
    assert(is_sqrtable());
    if (modulus() == 2 || value == 0) return *this;
    assert(prime::is_prime(modulus()));
    if (modulus() % 4 == 3) return pow((modulus() + 1) / 4);
    // get w=a^2-*this that is not sqrtable
    zint a, w;
    while (1) {
      a = random::rand(0, modulus() - 1);
      w = a * a - *this;
      if (!w.is_sqrtable()) break;
    }
    // x+y*sqrt(w)
    struct complex {
      zint x, y;
    };
    auto multiply_complex = [i2 = w](complex a, complex b) {
      return complex{a.x * b.x + a.y * b.y * i2, a.x * b.y + a.y * b.x};
    };
    complex res = {1, 0}, t = {a, 1};
    for (int k = (modulus() + 1) / 2; k; k >>= 1) {
      if (k & 1) res = multiply_complex(res, t);
      t = multiply_complex(t, t);
    }
    assert(res.x * res.x == *this);
    return std::min(res.x, -res.x);
  }

  bool operator==(const zint& z) const { return value == z.value; }
  bool operator!=(const zint& z) const { return value != z.value; }
  bool operator<(const zint& z) const { return value < z.value; }
  bool operator>(const zint& z) const { return value > z.value; }

  std::string to_string() const { return std::to_string(value); }

  friend std::ostream& operator<<(std::ostream& os, zint z) { return os << z.to_string(); }

  friend std::istream& operator>>(std::istream& is, zint& z) {
    int64_t t;
    is >> t;
    z = zint(t);
    return is;
  };

 private:
  type value;
  template <typename U>
  static typename std::enable_if_t<std::is_integral<U>::value, type> normalize(const U& x) {
    type value;
    if (-modulus() <= x && x < modulus())
      value = (type)x;
    else
      value = (type)(x % modulus());
    if (value < 0) value += modulus();
    return value;
  }
};
}  // namespace modular

template <typename T>
struct binomial {
  static T fac(int n) {
    static std::vector<T> f = {1};
    while (f.size() <= n) f.push_back(f.back() * f.size());
    return f[n];
  }
  static T inv(int n) {
    static std::vector<T> inv = {1, 1};
    while (inv.size() <= n)
      inv.push_back(-inv[T::modulus() % inv.size()] * (T::modulus() / inv.size()));
    return inv[n];
  }
  static T finv(int n) {
    static std::vector<T> finv = {1};
    while (finv.size() <= n) finv.push_back(finv.back() * inv(finv.size()));
    return finv[n];
  }
  static T C(int n, int m) {
    if (n < 0 || m < 0 || m > n) return T(0);
    return fac(n) * finv(m) * finv(n - m);
  }
  static T A(int n, int m) {
    if (n < 0 || m < 0 || m > n) return T(0);
    return fac(n) * finv(n - m);
  }
};

}  // namespace numerical

namespace polynomial {

template <typename T>
class poly : public std::vector<T> {
 public:
  using std::vector<T>::vector;

  poly(std::string s) {
    for (int i = 0; i < s.size();) {
      auto scan_num = [&]() -> long long {
        int sgn = 1;
        if (s[i] == '-') sgn = -1, ++i;
        if (s[i] == '+') sgn = 1, ++i;
        if (i == s.size() || !std::isdigit(s[i])) return sgn;
        long long num = 0;
        while (i < s.size() && std::isdigit(s[i])) num = num * 10 + s[i++] - '0';
        return sgn * num;
      };
      auto add_item = [&](size_t exponent, T coeff) {
        if (exponent >= this->size()) this->resize(exponent + 1);
        this->at(exponent) = coeff;
      };
      T coeff = scan_num();
      if (i == s.size() || s[i] != 'x')
        add_item(0, coeff);
      else {
        size_t exponent = 1;
        if (s[++i] == '^') {
          ++i;
          exponent = scan_num();
        }
        add_item(exponent, coeff);
      }
    }
  }

  int deg() const { return this->size() - 1; }
  poly operator-() const {
    poly ans = *this;
    for (auto& x : ans) x = -x;
    return ans;
  }
  T operator()(const T& x) const {
    T ans = 0;
    for (int i = this->size() - 1; i >= 0; --i) ans = ans * x + this->at(i);
    return ans;
  }
  T operator[](int idx) const {
    if (0 <= idx && idx < this->size()) return this->at(idx);
    return 0;
  }
  T& operator[](int idx) {
    if (idx >= this->size()) this->resize(idx + 1);
    return this->at(idx);
  }

  poly rev() const {
    poly res(*this);
    std::reverse(res.begin(), res.end());
    return res;
  }

  poly mulxk(size_t k) const {
    poly res = *this;
    res.insert(res.begin(), k, 0);
    return res;
  }
  poly divxk(size_t k) const {
    if (this->size() <= k) return {};
    return poly(this->begin() + k, this->end());
  }
  poly modxk(size_t k) const {
    k = std::min(k, this->size());
    return poly(this->begin(), this->begin() + k);
  }

  poly& operator*=(poly p);
  poly& operator+=(const poly& p) {
    this->resize(std::max(this->size(), p.size()));
    for (int i = 0; i < this->size(); ++i) this->at(i) += p[i];
    return this->normalize();
  }
  poly& operator-=(const poly& p) {
    this->resize(std::max(this->size(), p.size()));
    for (int i = 0; i < this->size(); ++i) this->at(i) -= p[i];
    return this->normalize();
  }
  poly& operator/=(const poly& p) {
    if (this->size() < p.size()) return *this = {};
    int len = this->size() - p.size() + 1;
    return *this = (this->rev().modxk(len) * p.rev().inv(len)).modxk(len).rev().normalize();
  }
  poly& operator%=(const poly& p) { return *this = (*this - (*this / p) * p).normalize(); }
  poly& operator*=(const T& x) {
    for (int i = 0; i < this->size(); ++i) this->at(i) *= x;
    return *this;
  }
  poly& operator/=(const T& x) { return *this *= x.inv(); }
  poly operator*(const poly& p) const { return poly(*this) *= p; }
  poly operator+(const poly& p) const { return poly(*this) += p; }
  poly operator-(const poly& p) const { return poly(*this) -= p; }
  poly operator/(const poly& p) const { return poly(*this) /= p; }
  poly operator%(const poly& p) const { return poly(*this) %= p; }
  poly operator*(const T& x) const { return poly(*this) *= x; }
  poly operator/(const T& x) const { return poly(*this) /= x; }

  // (quotient, remainder)
  std::pair<poly, poly> divmod(const poly& p) const {
    poly d = *this / p;
    return std::make_pair(d, (*this - d * p).normalize());
  }

  poly deriv() const {
    if (this->empty()) return {};
    poly res(this->size() - 1);
    for (int i = 0; i < this->size() - 1; ++i) {
      res[i] = this->at(i + 1) * (i + 1);
    }
    return res;
  }
  poly integr(T c = 0) const {
    poly res(this->size() + 1);
    for (int i = 0; i < this->size(); ++i) {
      res[i + 1] = this->at(i) / (i + 1);
    }
    res[0] = c;
    return res;
  }

  // mod x^k
  poly inv(int k = -1) const {
    if (!~k) k = this->size();
    poly res = {this->front().inv()};
    for (int len = 2; len < k * 2; len <<= 1) {
      res = (res * (poly{2} - this->modxk(len) * res)).modxk(len);
    }
    return res.modxk(k);
  }

  // mod x^k
  poly sqrt(int k = -1) const {
    if (!~k) k = this->size();
    poly res = {this->at(0).sqrt()};
    for (int len = 2; len < k * 2; len <<= 1) {
      res = (res + (this->modxk(len) * res.inv(len)).modxk(len)) / 2;
    }
    return res.modxk(k);
  }

  // mod x^k, a0=1 should hold
  poly log(int k = -1) const {
    assert(this->at(0) == 1);
    if (!~k) k = this->size();
    return (this->deriv() * this->inv(k)).integr().modxk(k);
  }
  // mod x^k, a0=0 should hold
  poly exp(int k = -1) const {
    assert(this->at(0) == 0);
    if (!~k) k = this->size();
    poly res = {1};
    for (int len = 2; len < k * 2; len <<= 1) {
      res = (res * (poly{1} - res.log(len) + this->modxk(len))).modxk(len);
    }
    return res.modxk(k);
  }
  // p^c mod x^k
  poly pow(int c, int k = -1) const {
    if (!~k) k = this->size();
    int i = 0;
    while (i < this->size() && !this->at(i)) ++i;
    if (i == this->size() || 1LL * i * c >= k) return {};
    T ai = this->at(i);
    poly f = this->divxk(i) * ai.inv();
    return (f.log(k - i * c) * c).exp(k - i * c).mulxk(i * c) * ai.pow(c);
  }

  // evaluate and interpolate
  struct product_tree {
    int l, r;
    std::unique_ptr<product_tree> lson = nullptr, rson = nullptr;
    poly product;
    product_tree(int l, int r) : l(l), r(r) {}

    static std::unique_ptr<product_tree> build(const std::vector<T>& xs,
                                               std::function<poly(T)> get_poly) {
      std::function<std::unique_ptr<product_tree>(int, int)> build = [&](int l, int r) {
        auto rt = std::make_unique<product_tree>(l, r);
        if (l == r) {
          rt->product = get_poly(xs[l]);
        } else {
          int mid = (l + r) >> 1;
          rt->lson = build(l, mid);
          rt->rson = build(mid + 1, r);
          rt->product = rt->lson->product * rt->rson->product;
        }
        return rt;
      };
      return build(0, xs.size() - 1);
    }
  };

  poly mulT(poly p) const {
    if (p.empty()) return {};
    return ((*this) * p.rev()).divxk(p.size() - 1);
  }

  std::vector<T> evaluate(std::vector<T> xs) const {
    if (this->empty()) return std::vector<T>(xs.size());
    std::unique_ptr<product_tree> rt = product_tree::build(xs, [&](T x) { return poly{1, -x}; });
    return evaluate_internal(xs, rt);
  }

  static poly interpolate(std::vector<T> xs, std::vector<T> ys) {
    assert(xs.size() == ys.size());
    if (xs.empty()) return {};
    std::unique_ptr<product_tree> rt = product_tree::build(xs, [&](T x) { return poly{1, -x}; });
    std::vector<T> coef = rt->product.rev().deriv().evaluate_internal(xs, rt);
    for (int i = 0; i < ys.size(); ++i) coef[i] = ys[i] * coef[i].inv();
    std::function<poly(product_tree*)> solve = [&](product_tree* rt) {
      if (rt->l == rt->r) {
        return poly{coef[rt->l]};
      } else {
        return solve(rt->lson.get()) * rt->rson->product.rev() +
               solve(rt->rson.get()) * rt->lson->product.rev();
      }
    };
    return solve(rt.get());
  }

  // a0=0 must hold
  poly cos(int k = -1) const {
    assert(this->at(0) == 0);
    if (!~k) k = this->size();
    T i = T::root().pow((T::modulus() - 1) / 4);
    poly x = *this * i;
    return (x.exp(k) + (-x).exp(k)) / 2;
  }

  // a0=0 must hold
  poly sin(int k = -1) const {
    assert(this->at(0) == 0);
    if (!~k) k = this->size();
    T i = T::root().pow((T::modulus() - 1) / 4);
    poly x = *this * i;
    return (x.exp(k) - (-x).exp(k)) / (i * 2);
  }

  // a0=0 must hold
  poly tan(int k = -1) const { return this->sin(k) / this->cos(k); }

  poly acos(int k = -1) const {
    const poly& x = *this;
    return (-x.deriv() * (poly{1} - x * x).sqrt().inv()).integr();
  };
  poly asin(int k = -1) const {
    const poly& x = *this;
    return (x.deriv() * (poly{1} - x * x).sqrt().inv()).integr();
  };
  poly atan(int k = -1) const {
    const poly& x = *this;
    return (x.deriv() * (poly{1} + x * x).inv()).integr();
  };

  friend std::ostream& operator<<(std::ostream& os, poly p) {
    os << "{";
    for (auto x : p) os << x() << " ";
    os << "}";
    return os;
  }

 private:
  poly& normalize() {
    while (this->size() && !this->back()) this->pop_back();
    return *this;
  }

  std::vector<T> evaluate_internal(std::vector<T>& xs, std::unique_ptr<product_tree>& rt) const {
    std::vector<T> res(xs.size());
    xs.resize(std::max(xs.size(), this->size()));
    std::function<void(product_tree*, poly)> solve = [&](product_tree* rt, poly p) {
      p = p.modxk(rt->r - rt->l + 1);
      if (rt->l == rt->r) {
        if (rt->l < res.size()) res[rt->l] = p.front();
      } else {
        solve(rt->lson.get(), p.mulT(rt->rson->product));
        solve(rt->rson.get(), p.mulT(rt->lson->product));
      }
    };
    solve(rt.get(), this->mulT(rt->product.inv(xs.size())));
    return res;
  }
};

template <typename T>
class dft {
 public:
  static const bool use_fast_trans = true;
  static void trans(poly<T>& p) {
    assert(__builtin_popcount(p.size()) == 1);
    if constexpr (use_fast_trans) {
      dif(p);
    } else {
      bit_reverse(p);
      dit(p);
    }
  }

  static void inv_trans(poly<T>& p) {
    assert(__builtin_popcount(p.size()) == 1);
    if constexpr (use_fast_trans) {
      dit(p);
    } else {
      trans(p);
    }
    reverse(p.begin() + 1, p.end());
    T inv = T(p.size()).inv();
    for (T& x : p) x *= inv;
  }

  // should call dit after dif
  static void dit(poly<T>& p) {
    for (int len = 1; len < p.size(); len <<= 1) {
      auto sub_w = get_subw(len * 2);
      for (auto sub_p = p.begin(); sub_p != p.end(); sub_p += 2 * len)
        for (int i = 0; i < len; ++i) {
          T u = sub_p[i], v = sub_p[i + len] * sub_w[i];
          sub_p[i] = u + v;
          sub_p[i + len] = u - v;
        }
    }
  }

  static void dif(poly<T>& p) {
    for (int len = p.size() / 2; len >= 1; len >>= 1) {
      auto sub_w = get_subw(len * 2);
      for (auto sub_p = p.begin(); sub_p != p.end(); sub_p += 2 * len)
        for (int i = 0; i < len; ++i) {
          T _sub_pi = sub_p[i];
          sub_p[i] += sub_p[i + len];
          sub_p[i + len] = (_sub_pi - sub_p[i + len]) * sub_w[i];
        }
    }
  }

 private:
  typename std::vector<T>::iterator static get_subw(int len) {
    static std::vector<T> w = {0, 1};
    while (w.size() <= len) {
      T e[] = {1, T::root().pow((T::modulus() - 1) / w.size())};
      w.resize(w.size() * 2);
      for (int i = w.size() / 2; i < w.size(); ++i) w[i] = w[i / 2] * e[i & 1];
    }
    return w.begin() + len;
  }
};

template <typename T>
poly<T>& poly<T>::operator*=(poly<T> p) {
  if (this->empty() || p.empty()) return *this = {};
  constexpr int small_size = 128;
  if (this->size() < small_size || p.size() < small_size) {
    poly<T> t(this->size() + p.size() - 1);
    for (int i = 0; i < this->size(); i++)
      for (int j = 0; j < p.size(); j++) t[i + j] += this->at(i) * p[j];
    return *this = t;
  }
  int len = 1 << (std::__lg(this->deg() + p.deg()) + 1);
  this->resize(len);
  p.resize(len);
  dft<T>::trans(*this);
  dft<T>::trans(p);
  for (int i = 0; i < len; ++i) this->at(i) *= p[i];
  dft<T>::inv_trans(*this);
  return this->normalize();
}

namespace utils {

// calc [x^n]g(x) than f(g(x))=h(x)
// f[0]=0 should hold
template <typename T>
T lagrange_inversion(poly<T> f, int n, poly<T> h = poly<T>{0, 1}) {
  assert(f[0] == 0);
  return (f.divxk(1).pow(n).inv() * h.deriv())[n - 1] * zayin::numerical::binomial<T>::inv(n);
}

}  // namespace utils

}  // namespace polynomial

namespace graph_theory {

template <typename T>
struct weighted_edge {
  int to;
  T weight;
  weighted_edge(int to, T weight) : to(to), weight(weight) {}
};

template <typename T>
struct graph {
  const int n;
  std::vector<T> edges;
  std::vector<std::vector<int>> g;
  graph(int n) : n(n), g(n) {}

  void add_edge(int u, T edge) {
    assert(0 <= u && u < n);
    g[u].push_back(edges.size());
    edges.push_back(edge);
  }
};

namespace flow {

template <typename T, typename std::enable_if_t<std::is_arithmetic<T>::value, int> = 0>
struct maxflow : public graph<weighted_edge<T>> {
  std::vector<int> cur, h;
  maxflow(int n) : graph<weighted_edge<T>>(n), cur(n), h(n) {}

  void add_edge(int u, int v, T w) {
    graph<weighted_edge<T>>::add_edge(u, weighted_edge<T>(v, w));
    graph<weighted_edge<T>>::add_edge(v, weighted_edge<T>(u, 0));
  }

  bool bfs(int s, int t) {
    h.assign(this->n, -1);
    std::queue<int> q;
    h[s] = 0;
    q.push(s);
    while (!q.empty()) {
      int u = q.front();
      q.pop();
      for (int i : this->g[u]) {
        auto [v, w] = this->edges[i];
        if (w > 0 && !~h[v]) {
          h[v] = h[u] + 1;
          if (v == t) return true;
          q.push(v);
        }
      }
    }
    return false;
  }

  T dfs(int u, int t, T f) {
    if (u == t) return f;
    T r = f;
    for (int& i = cur[u]; i < this->g[u].size(); ++i) {
      int j = this->g[u][i];
      auto [v, w] = this->edges[j];
      if (w > 0 && h[v] == h[u] + 1) {
        auto a = dfs(v, t, std::min(r, w));
        this->edges[j].weight -= a;
        this->edges[j ^ 1].weight += a;
        r -= a;
        if (!r) break;
      }
    }
    return f - r;
  }

  T max_flow(int s, int t) {
    T res = 0;
    while (bfs(s, t)) {
      cur.assign(this->n, 0);
      res += dfs(s, t, std::numeric_limits<T>::max());
    }
    return res;
  }
};

}  // namespace flow

namespace shortest_path {

template <typename T, typename std::enable_if_t<std::is_arithmetic<T>::value, int> = 0>
struct dijkstra : public graph<weighted_edge<T>> {
  static constexpr T inf = std::numeric_limits<T>::max();

  std::vector<T> dis, reachable;

  dijkstra(int n) : graph<weighted_edge<T>>(n), dis(n), reachable(n) {}

  void init(int vs) {
    struct item {
      int u;
      T d;
      bool operator<(const item& I) const { return d > I.d; }
    };
    dis.assign(this->n, inf);
    reachable.assign(this->n, false);
    std::priority_queue<item> q;
    q.push({vs, dis[vs] = 0});
    while (q.size()) {
      auto [u, d] = q.top();
      q.pop();
      if (dis[u] < d) continue;
      reachable[u] = true;
      for (int i : this->g[u]) {
        auto [v, w] = this->edges[i];
        if (dis[v] <= dis[u] + w) continue;
        q.push({v, dis[v] = dis[u] + w});
      }
    }
  }
};

}  // namespace shortest_path

}  // namespace graph_theory

namespace data_structure {
namespace segment_tree {

template <typename T>
class segment_tree {
 public:
  segment_tree(int n_) {
    for (n = 1; n < n_; n <<= 1) {
    }
    tree.resize(n << 1);
  }
  template <typename U>
  segment_tree(
      std::vector<U> a, std::function<T(U)> trans = [](const U& u) { return T(u); }) {
    for (n = 1; n < (int)a.size(); n <<= 1) {
    }
    tree.resize(n << 1);
    for (int i = 0; i < a.size(); ++i) tree[n + i] = trans(a[i]);
    for (int i = n; --i;) tree[i] = T::merge(tree[i << 1], tree[i << 1 | 1]);
  }

  void change(int i, const T& t) {
    assert(0 <= i && i < n);
    tree[i += n] = t;
    while (i >>= 1) tree[i] = T::merge(tree[i << 1], tree[i << 1 | 1]);
  }

  // [l,r]
  T composite(int l, int r) {
    assert(0 <= l && l < n);
    assert(0 <= r && r < n);
    if (l > r) return T();
    T prodL, prodR;
    for (l += n, r += n + 1; l < r; l >>= 1, r >>= 1) {
      if (l & 1) prodL = T::merge(prodL, tree[l++]);
      if (r & 1) prodR = T::merge(tree[--r], prodR);
    }
    return T::merge(prodL, prodR);
  }

 private:
  int n;
  std::vector<T> tree;
};

class lazy_segment_tree;

};  // namespace segment_tree

namespace lichao_tree {

constexpr long long inf32 = 1e9;
constexpr long long inf64 = 1e18;

struct line {
  long long k, b;
  line(long long k = inf32, long long b = inf64) : k(k), b(b) {}
  long long operator()(long long x) const { return k * x + b; }
  friend std::ostream& operator<<(std::ostream& os, line line) {
    return os << "(" << line.k << "," << line.b << ")";
  }
  std::string to_string() const {
    std::stringstream os;
    os << *this;
    return os.str();
  }
};

// T is a monotonic function, and two different T should have most one intersection
template <typename T, typename cmp = std::less<decltype(T().operator()(0))>>
class lichao_tree {
 public:
  using type = decltype(T().operator()(0));
  static_assert(std::is_integral<type>::value);
  static const type type_min =
      std::min(std::numeric_limits<type>::min(), std::numeric_limits<type>::max(), cmp());
  static const type type_max =
      std::max(std::numeric_limits<type>::min(), std::numeric_limits<type>::max(), cmp());

  lichao_tree(type domain_left, type domain_right, size_t maxsize = 0)
      : domain_left(domain_left), domain_right(domain_right) {
    assert(domain_left < domain_right);
    rt = -1;
    best.reserve(maxsize);
    ls.reserve(maxsize);
    rs.reserve(maxsize);
  }

  void insert(T f) { rt = insert(rt, domain_left, domain_right, f); }
  type query(type x) { return query(rt, domain_left, domain_right, x); }

 private:
  int newnode(T f) {
    int id = best.size();
    best.emplace_back(f);
    ls.push_back(-1);
    rs.push_back(-1);
    return id;
  }
  int insert(int rt, type l, type r, T f) {
    if (!~rt) {
      rt = newnode(f);
    } else {
      type mid = l + (r - l) / 2;
      if (cmp()(f(mid), best[rt](mid))) std::swap(best[rt], f);
      if (cmp()(f(l), best[rt](l))) ls[rt] = insert(ls[rt], l, mid, f);
      if (cmp()(f(r), best[rt](r))) rs[rt] = insert(rs[rt], mid + 1, r, f);
    }
    return rt;
  }

  type query(int rt, type l, type r, type x) {
    if (!~rt) return type_max;
    if (l == r) return best[rt](x);
    type mid = l + (r - l) / 2;
    if (x <= mid)
      return std::min(best[rt](x), query(ls[rt], l, mid, x), cmp());
    else
      return std::min(best[rt](x), query(rs[rt], mid + 1, r, x), cmp());
  }

  int rt;
  std::vector<T> best;
  std::vector<int> ls, rs;
  const type domain_left;
  const type domain_right;
};

}  // namespace lichao_tree

namespace linear_basis {
template <typename T, int width = sizeof(T) * 8>
class linear_basis {
 public:
  static_assert(std::is_integral<T>::value);
  linear_basis() { bases_.resize(width); }
  size_t basic_count() { return basic_count_; }
  std::vector<T> bases() { return bases_; }
  bool insert(T x) {
    for (int i = width - 1; i >= 0; --i) {
      if (~x >> i & 1) continue;
      if (bases_[i]) {
        x ^= bases_[i];
      } else {
        ++basic_count_;
        bases_[i] = x;
        return true;
      }
    }
    return false;
  }

  T maximum() const {
    T res = 0;
    for (int i = width - 1; i >= 0; --i) res = std::max(res, res ^ bases_[i]);
    return res;
  }

  friend linear_basis operator+(linear_basis lhs, linear_basis rhs) {
    for (T x : lhs.bases()) rhs.insert(x);
    return rhs;
  }

 private:
  size_t basic_count_ = 0;
  std::vector<T> bases_;
};

}  // namespace linear_basis

namespace string {

template <size_t alphabet_size = 26>
class suffix_automaton {
 public:
  struct node {
    std::array<int, alphabet_size> to;
    int link, len, count;

    explicit node(int len = 0) : len(len), link(-1), count(0) { to.fill(-1); }
    explicit node(int len, int link, const std::array<int, alphabet_size>& to)
        : len(len), link(link), to(to), count(0) {}
  };

  suffix_automaton() : nodes() { newnode(); }

  explicit suffix_automaton(const std::string& s) : suffix_automaton() { insert(s); }

  void insert(const std::string& s) {
    nodes.reserve(size() + s.size() * 2);
    int last = 0;
    for (int i = 0; i < s.size(); ++i) {
      last = extend(last, s[i] - 'a');
    }
  }

  int extend(int k, int c) {
    if (~nodes[k].to[c] && nodes[nodes[k].to[c]].len == nodes[k].len + 1) {
      return nodes[k].to[c];
    }
    int leaf = newnode(nodes[k].len + 1);
    for (; ~k && !~nodes[k].to[c]; k = nodes[k].link) nodes[k].to[c] = leaf;
    if (!~k) {
      nodes[leaf].link = 0;
    } else {
      int p = nodes[k].to[c];
      if (nodes[k].len + 1 == nodes[p].len) {
        nodes[leaf].link = p;
      } else {
        int np = newnode(nodes[k].len + 1, nodes[p].link, nodes[p].to);
        nodes[p].link = nodes[leaf].link = np;
        for (; ~k && nodes[k].to[c] == p; k = nodes[k].link) nodes[k].to[c] = np;
      }
    }
    return leaf;
  }

  void mark_count(const std::string& s) {
    int k = 0;
    for (char c : s) {
      k = nodes[k].to[c - 'a'];
      assert(~k);
      ++nodes[k].count;
    }
  }

  void addup_count() {
    std::vector<int> ids(size()), bucket(size());
    for (int i = 0; i < size(); ++i) ++bucket[nodes[i].len];
    for (int i = 1; i < bucket.size(); ++i) bucket[i] += bucket[i - 1];
    for (int i = 0; i < size(); ++i) ids[--bucket[nodes[i].len]] = i;
    for (int i = size() - 1; i; --i) nodes[nodes[ids[i]].link].count += nodes[ids[i]].count;
  }

  const node& operator[](int v) const { return nodes[v]; }

  int maxlen(int v) const { return nodes[v].len; }

  int minlen(int v) const { return v ? nodes[nodes[v].link].len + 1 : 0; }

  int size() const { return nodes.size(); }

  std::string to_string() const {
    std::ostringstream os;
    std::function<void(int, std::string)> travel = [&](int k, std::string s) {
      if (!~k) return;
      os << k << ": " << s << " ~ " << nodes[k].count << "\n";
      for (int c = 0; c < alphabet_size; ++c) travel(nodes[k].to[c], s + (char)(c + 'a'));
    };
    travel(0, "");
    return os.str();
  }

 private:
  std::vector<node> nodes;

  template <typename... Args>
  int newnode(Args... args) {
    int res = nodes.size();
    nodes.push_back(node{args...});
    return res;
  }
};
}  // namespace string

}  // namespace data_structure

}  // namespace zayin

namespace io = zayin::io;
using zint = zayin::numerical::modular::zint<zayin::numerical::modular::modulus_type>;
namespace polynomial = zayin::polynomial;
using poly = polynomial::poly<zint>;
using binomial = zayin::numerical::binomial<zint>;
namespace prime = zayin::numerical::prime;

#pragma endregion

int main() {
  std::ios::sync_with_stdio(false);
  std::cin.tie(0);
  std::cout.tie(0);
  std::string s;
  std::cin >> s;
  int n = s.size();
  auto solve1 = [&]() {
    int ans = 0;
    for (int i = 1; i < n; ++i) {
      std::string a = s.substr(0, i), b = s.substr(i);
      std::vector<std::vector<int>> dp(a.size() + 1, std::vector<int>(b.size() + 1));
      for (int i = 1; i <= a.size(); ++i)
        for (int j = 1; j <= b.size(); ++j) {
          if (a[i - 1] == b[j - 1]) {
            dp[i][j] = dp[i - 1][j - 1] + 1;
          } else {
            dp[i][j] = std::max(dp[i - 1][j], dp[i][j - 1]);
          }
        }
      ans = std::max(ans, dp[a.size()][b.size()] * 2);
    }
    return ans;
  };
  auto solve2 = [&]() {
    int ans = 0;
    for (int i = 1; i < n; ++i) {
      for (int j = i + 1; j < n; ++j) {
        std::string a = s.substr(0, i), b = s.substr(i, j - i), c = s.substr(j);
        std::vector<std::vector<std::vector<int>>> dp(
            a.size() + 1,
            std::vector<std::vector<int>>(b.size() + 1, std::vector<int>(c.size() + 1)));
        for (int i = 1; i <= a.size(); ++i)
          for (int j = 1; j <= b.size(); ++j) {
            for (int k = 1; k <= c.size(); ++k) {
              if (a[i - 1] == b[j - 1] && b[j - 1] == c[k - 1]) {
                dp[i][j][k] = dp[i - 1][j - 1][k - 1] + 1;
              } else {
                dp[i][j][k] = std::max({dp[i - 1][j][k], dp[i][j - 1][k], dp[i][j][k - 1]});
              }
            }
          }
        ans = std::max(ans, dp[a.size()][b.size()][c.size()] * 3);
      }
    }
    return ans;
  };

  auto solve3 = [&]() {
    int ans = 0;
    constexpr int bench_len = 16;
    for (int l = 0; l < n; l += bench_len) {
      int r = std::min(l + bench_len, n);
      for (int mask = 1; mask < (1 << (r - l)); ++mask) {
        std::string t;
        for (int i = l; i < r; ++i)
          if (mask >> (i - l) & 1) t.push_back(s[i]);
        int len = 0, j = 0;
        for (int i = 0; i < n; ++i) {
          if (s[i] == t[j]) {
            ++len;
            j = (j + 1) % t.size();
          }
        }
        int cnt = len / t.size();
        if (cnt > 1) ans = std::max<int>(ans, cnt * t.size());
      }
    }
    return ans;
  };

  std::cout << std::max({solve1(), solve2(), solve3()}) << std::endl;
  return 0;
}